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We study in detail the active Ising model, a stochastic lattice gas where collective motion emerges 
from the spontaneous breaking of a discrete symmetry. On a 2d lattice, active particles undergo a 
diffusion biased in one of two possible directions (left and right) and align ferromagnetically their 
direction of motion, hence yielding a minimal flocking model with discrete rotational symmetry. 
We show that the transition to collective motion amounts in this model to a bona fide liquid-gas 
phase transition in the canonical ensemble. The phase diagram in the density/velocity param¬ 
eter plane has a critical point at zero velocity which belongs to the Ising universality class. In 
the density/temperature ‘canonical’ ensemble, the usual critical point of the equilibrium liquid-gas 
transition is sent to infinite density because the different symmetries between liquid and gas phases 
preclude a supercritical region. We build a continuum theory which reproduces qualitatively the 
behavior of the microscopic model. In particular we predict analytically the shapes of the phase 
diagrams in the vicinity of the critical points, the binodal and spinodal densities at coexistence, and 
the speeds and shapes of the phase-separated profiles. 
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I. INTRODUCTION 

Active matter systems, defined as large assemblies of 
interacting particles consuming energy to self-propel, ex¬ 
hibit a variety of elaborate collective behaviors. Among 
them, collective motion—a term referring to the coherent 
displacement of large groups of individuals over length 
scales much larger than their individual size—has played 
a leading role in active matter studies. It can be ob¬ 
served in a wide range of biological systems such as bird 
flocks [1], fish schools m 0, bacterial swarms 0 0, 
actin [6| or microtubule [7] motility assays, but also in 
inert matter that is artificially self-propelled, for exam¬ 
ple in assemblies of vibrated polar disks [8], rolling col¬ 
loids [9] or self-propelled liquid droplets m- 

On the theoretical side, the transition to collec¬ 
tive motion—hereafter referred to as the “flocking” 
transition—has attracted the attention of the physics 
community because simple models have proved useful to 
describe its generic properties, highlighting the possibil¬ 
ity of universal behaviours. The model introduced by 
Vicsek and collaborators two decades ago m is proto¬ 
typical of this line of research, containing only two ingre¬ 
dients: self-propulsion at a constant speed and aligning 
interactions. It has often been described as a dynamical 
XY model [TO] since the alignment of the particle direc¬ 
tions of motion resemble the ferromagnetic alignment of 
XY spins. 

The phenomenology of the Vicsek model is now well es- 
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tablished mmsun]. When decreasing the noise on the 
aligning interaction, or increasing the density, a transi¬ 
tion takes place from a disordered gas into an ordered 
state of collective motion. Between these two homoge¬ 
neous phases lays a region of parameter space where par¬ 
ticles gather in dense ordered bands travelling in a dilute 
disordered background. These bands, which are a robust 
feature of flocking models [iniiisHsni, are a signature of 
the first-order nature of the transition, together with in- 
termittency, metastability and hysteresis [131 [H]. Unfor¬ 
tunately, they are seen only in large systems and strong 
finite size effects render the numerical study of the Vicsek 
model (VM) very costly in computing power. 

To overcome these numerical difficulties and gain 
more insight into the flocking transition, a number of 
analytical approaches have been followed. Hydrody¬ 
namic equations for Vicsek-like models have been ei¬ 
ther derived by coarse-graining [Huni or proposed phe¬ 
nomenologically [121 IISI 121] • These equations predict 
phase diagrams in qualitative agreement with the micro¬ 
scopic models, including the existence of inhomogeneous 
bands [HUa [111120]. Their analytical study is however 
so complicated that little can be done beyond working 
with their linearized version. Nevertheless, some progress 
was made to account for the long range order and gi¬ 
ant density fluctuations observed in the ordered phase of 
the Vicsek model m- Interestingly, it was also recently 
shown that all hydrodynamic equations derived for polar 
flocking models [m [1111211122] admit the same family 
of Id propagative solutions [23]. A complete analytical 
study of the Vicsek model, from micro to macro, however 
remains elusive. 

An alternative strategy to gain insight into the flock¬ 
ing transition relied on the introduction of an Active Ising 
Model (AIM) [22] which circumvents both the numerical 
and analytical pitfalls of the Vicsek model. Using non- 
equilibrium versions of ferromagnetic models has indeed 
often proven a useful strategy [TT] |231[26] . The AIM, 
which we study in detail in this paper, contains the two 
key ingredients for flocking: self-propulsion and aligning 
interactions. The continuous rotational symmetry of the 
Vicsek model is however replaced by a discrete symme¬ 
try; In the AIM, particles diffuse in the 2d plane but 
are self-propelled in only one of two possible directions 
(left or right). It is thus akin to a dynamical Ising model 
where particles have a discrete rotational symmetry. The 
AIM is found to have a simpler, more tractable, behavior 
than the Vicsek-like models with continuous symmetry 
while still retaining a large part of their physics. Using 
a lattice-based model also simplifies both numerical and 
analytical studies. 

After introducing the model in section [Til w e present a 
numerical study of the 2d AIM in section ]lTI| Our main 
conclusion is that the transition in the AIM amounts 
to a liquid-gas transition in the canonical ensemble. At 
fixed orientational noise, the system can be in two “pure” 
states: a disordered gas or an ordered liquid, the latter 
leading to a collective migration of all particles to the 
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FIG. 1. Sketch of the two possible actions and their rates of 
occurrence. The ferromagnetic interaction between particles 
is purely on-site and particles diffuse freely. Beyond the biased 
diffusion shown here, particles also hop symmetrically up or 
down, with equal rates D in both directions. 


left or to the right. When constraining the system’s den¬ 
sity to he between two ’spinodal lines’, no homogeneous 
phase can be observed and the system phase-separates, 
with an ordered travelling liquid band coexisting with a 
disordered gas background. A key difference with the 
usual equilibrium liquid-gas transition is that liquid and 
gas have different symmetries; A supercritical region is 
thus prohibited since one has to break a symmetry to take 
the system from a gas to a liquid state, which explains 
the atypical shape of the phase diagram. 

In section [TV| we complement our numerical approach 
by deriving a set of hydrodynamic equations for the dy¬ 
namics of the local density and magnetisation fields. In¬ 
terestingly, a simple mean-field theory wrongly predicts 
a continuous transition, failing to account for the phase- 
separated profiles. A refined mean-field model, taking 
into account the fluctuations of the density and magneti¬ 
sation fields, reproduces qualitatively the phenomenol¬ 
ogy of the AIM. In section |V| we use the hydrodynamic 
equations to compute at large densities the shape of the 
phase-separated profiles, the coexisting densities, the ve¬ 
locity of the liquid domain and account for the finite-size 
scalings observed in the microscopic model. Finally, we 
argue in section [VT] in favor of the robustness of our re¬ 
sults by considering an off-lattice version of the model 
and different boundary conditions. 


II. DEFINITION OF THE MODEL 


We consider N particles moving on a 2D lattice of 
Lx X Ly sites with periodic boundary conditions. Each 
particle carries a spin ±I and there are no excluded vol¬ 
ume interactions between the particles: there can thus 
be an arbitrary number nf^ of particles with spins ±I on 
each site i = (UG 2 )- The local densities and magnetiza¬ 
tions are then defined as pi = nj' and rrii = n'l —n^. 
We consider a continuous-time Markov process in which 
particles can both flip their spins and hop to neighbor¬ 
ing lattice sites at rates that depend on their spins (see 
Fig. [^. The hopping and flipping rates, detailed in the 
next subsections, are such that our model is endowed 
with self-propulsion and inter-particle alignment, hence 
consituting a flocking model with discrete symmetry. 
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FIG. 2. A loop of four configurations involving 2 particles on 
2 sites breaking Kolmogorov’s criterion [^7] showing that the 
system does not satisfy detailed balance even when e = 0. 
The numbers associated to the arrows are the transition rates 
for £ = 0. The product of the transition rates along Ci 
C 2 ^ C 3 ^ C 4 ^ Cl (left to right) is , whereas the 

reverse order (right to left) yields 2D^e~^. 


A. Alignment: Fully connected Ising models 
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A particle with spin s on site i flips its spin at rate 


hF(s ^ —s) = yexp s/3—^ , ( 1 ) 

where (3 = 1/T plays the role of an inverse temperature. 
These rates satisfy detailed balance with respect to an 
equilibrium distribution P oc ex.p[—/3H] where H is the 
sum over the L^Ly lattice sites of the Hamiltonians of 
fully connected Ising models: 


sites i j=l kj^j 


sites 1 


mt 


, [ 2 ^ 


( 2 ) 


The first sum runs over the lattice site index i = ( 21 ^ 2 ), 
the next two over the particles j, k present on site i, and 
Sj = ±lis the value of spin j. (The factor 1/2 simply 
avoids double counting.) The rate 7 can always be ab¬ 
sorbed in a change of time unit so that we take 7 = 1 , 
silently omitting it from now on. 

This interaction is purely local: particles only align 
with other particles on the same site and, without parti¬ 
cle hopping, the model amounts to independent fully 
connected Ising models. The factor 1/pi in W makes the 
Hamiltonian H extensive with N and keeps the inter¬ 
action rates bounded: the rate W{s ^ —s) at which a 
particle of spin s flips its spin varies between exp(—/3) if 
all the other particles on the same site have spins s to 
exp[/3(l — 2/Pi)] if they all have spins —s. 


B. Self-propulsion: Biased diffusion 

Particles also undergo free diffusion on the lattice, with 
a left/right bias depending on the sign of their spins: a 
particle with spin s hops with rate se) to its right, 

D (1 — se) to its left, and D in both the up and down 
directions. There is thus a mean drift, which plays the 
role of self-propulsion, with particles of spins ±1 moving 
along the horizontal axis with an average velocity P2De. 


FIG. 3. (a-c) Examples of density profiles (green upper line) 
and magnetization profiles (blue lower line) averaged along 
vertical direction for the three phases, (a) Disordered gas, 
/3 = 1.4, po = 2. (b) Polar liquid, ^ = 2, po = 7- (c) 

Liquid-gas coexistence, /3 = 1.6, po = 5. (d) 2d snapshot 
corresponding to (c). (for all figures D = 1, £ = 0.9) 


The model is designed to have the self-propulsion en¬ 
tering in a minimal and tunable way through the param¬ 
eter £. Importantly, the limit of vanishing self-propulsion 
£ ^ 0 is well-defined because the spins still diffuse on the 
lattice. This dynamics should thus allow us to interpole 
continuously between ‘totally self-propelled’ {e = 1 ), self- 
propelled {e g] 0, 1[), weakly self-propelled {e ^ 1/L) and 
purely diffusive {e = 0 ) particles. 


This differs from the Vicsek model where the zero- 
velocity limit corresponds to immobile particles under¬ 
going an equilibrium dynamics resembling that of the 
XY model, with a quenched disorder on the bonds (only 
particles closer than a fixed distance interact). 

Let us note, however, that even when 5 = 0 the model 
is not at equilibrium i.e. it does not satisfy detailed 
balance with respect to any distribution. This is easily 
shown using Kolmogorov’s criterion m- In Fig. we 
exhibit a loop of four configurations such that the prod¬ 
ucts of the transition rates for visiting the loop in one 
order, Ci ^ C 2 ^ C 3 ^ C 4 —> Ci, and the reverse order 
are different, whence a violation of detailed balance. To 
make the 5 = 0 limit an equilibrium dynamics, one strat¬ 
egy could be to choose hopping rates satisfying detailed 
balance with respect to the Hamiltonian H defined in Q , 
replacing D by D ex.-p{—/3AH/2)). The steady-state dis¬ 
tribution would however be factorized and not very in¬ 
teresting. An alternative would be to further add to (© 
nearest neighbours interactions but we have not followed 
this cumbersome path here. Actually, as we show in sec¬ 
tion HIB , this microscopic irreversibility when 5 = 0 is 
irrelevant at large scales and we recover in this limit a 
phase transition belonging to the Ising universality class. 
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FIG. 4. Phase-separated density (left) and magnetization 
(right) profiles as the density increases. Parameters: /3 = 2, 
D = 1, £ = 0.9, system of size 800x100. The profiles have 
been averaged over time and along the y axis. 


C. Simulations 




FIG. 5. Phase diagrams of the AIM. The red and blue lines de¬ 
limit the region of existence of phase-separated profiles. Left: 
Parameter spaces (T = 1/^, po) for D = 1. Red and blue co¬ 
existence lines correspond to £ = 0.9 while the green dashed 
line indicates the critical points at £ = 0. Right: Parameter 
space (s, Po) for D = 1, ^ = 1.9. At e = 0 we recover critical 
point in the Ising universality class. 


To simulate the dynamics of the model, we used a 
random-sequential-update algorithm. We discretized the 
time in small time-steps At. A particle is then chosen at 
random; it flips its spin s with probability W{s ^ —s) At, 
hops upwards or downwards with probabilities DAt, to 
its right or to its left neighboring sites with probabilities 
D{1 A se)At. Finally, it does nothing with probability 
1 — [4:D -\-W{s ^ —s)]At. Time is then incremented by 
At/N and we iterate up to some final time. In practice 
we used At = [AD ^ exp{/3)]~^ to minimize the probabil¬ 
ity that nothing happens while keeping all probabilities 
smaller than one. 

Note that this algorithm does not allow a particle to be 
updated twice (on average) during At and is thus an ap¬ 
proximation of our continuous-time Markov process. We 
also used continuous-time simulations, associating clocks 
to each particle or each site and pulling updating times 
from the corresponding exponential laws. In practice we 
did not find any difference in the simulation results but 
the continuous time simulations were often slower so that 
we mostly used the random sequential update algorithm. 

In most of this article we use simulation boxes with 
Lx X Ly lattice sites and periodic boundary conditions. In 
section [VI A| we discuss what happens for closed bound¬ 
ary conditions. 


III. A LIQUID-GAS PHASE TRANSITION 

We explored the phase diagram using three control pa¬ 
rameters: the temperature T = the average density 
po = N/{LxLy)^ and the self-propulsion ‘speed’ 5. Doing 
so, we observed three different phases shown in Fig. 
For 5 7 ^ 0, at high temperatures and low densities, the 
particles fail to organize and we observe a homogeneous 
gas of particles with local magnetization {rui) ~ 0. On 
the contrary, for large densities and small temperatures, 
the particles move collectively either to the right or to the 
left, forming a polar liquid state with (m^) = mo 7 ^ 0 . 
For intermediate densities, when po G [p^(T, e), p^(T, e)], 



FIG. 6. Schematic picture of the differences between the 
phase diagrams of the passive and active liquid gas transi¬ 
tion. In the active case, because the liquid and the gas have 
different symmetries, the critical point is sent to p = 00 , thus 
suppressing the supercritical region. 


the system phase separates into a band of polar liquid 
traveling to the left or to the right through a disordered 
gaseous background. 

The lines p^(T, e) and p^(T, e) both delimit the do¬ 
main of existence of the phase-separated profiles and play 
the role of coexistence lines: As shown in Fig. for 
all phase-separated profiles at fixed T, e, the densities 
in the gas and liquid part of the profiles are pg and p^, 
respectively. Correspondingly, the magnetization are 0 
and mi{T^e) 7 ^ 0. Thus, varying the density po at con¬ 
stant temperature and propulsion speed solely changes 
the width of the liquid band. Consequently, in the phase 
coexistence region, the lever rule can be used to deter¬ 
mine the liquid fraction 4> in the same way as for an 
equilibrium liquid-gas phase transition in the canonical 
ensemble: 


$ = ^^ 
Pi - Pg 


(3) 


As we shall see below, this analogy goes beyond the 
sole shape of the phase-separated profiles and the phase- 
transition to collective motion of the active Ising model 
is best described as a liquid-gas phase transition rather 
than an order-disorder one. 
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A. Temperature-density ‘canonical’ ensemble 

The phase diagram in the (T, po) parameter plane, 
computed for 5 = 0.9, is shown in the left panel of Fig.|^ 
While the general structure of the phase diagram, with 
a gas phase, a liquid phase, and a coexistence region, is 
reminiscent of an equilibrium liquid-gas phase diagram, 
the shapes of the transition lines are unusual. This dif¬ 
ference can be understood using a symmetry argument. 
Since the disordered gas and the polar liquid have differ¬ 
ent symmetries, the system cannot continuously trans¬ 
form from one homogeneous phase to the other without 
crossing a transition line. There is thus no super-critical 
region and the critical point is sent to Tc = 1 and pc = 00 . 
(See Fig. for a schematic picture.) 

This symmetry argument should be rather general for 
flocking transitions separating a disordered state and a 
symmetry-broken state of collective motion. Indeed, in 
Vicsek-like models, where the role of the inverse temper¬ 
ature is played by the noise intensity, the phase diagrams 
are qualitatively similar to the one shown in Fig.|^ This 
is true both for the full phase diagram recently computed 
in M as well as for earlier results m, for a slightly dif¬ 
ferent kinetic model and its hydrodynamic theory m, 
but also for an active nematic Vicsek-like model [28] and 
a hydrodynamic theory of self-propelled rods [29] . 


B. Velocity-density ensemble 


Conversely, one can change the strength of the self¬ 
propulsion 5 while keeping the temperature fixed. Again, 
one obtains a phase diagram with three regions. The dif¬ 
ference with the canonical ensemble is that in this pa¬ 
rameter plane, the two coexistence lines merge at 5 = 0, 
where self-propulsion vanishes, yielding a critical point at 
a finite density p*(T) (See the right panel of Fig. |§. The 
curve p*(T) is reported in the left panel of Fig. [ 5 ] and sat¬ 
isfies p*(T) G [pg{T,e), p£{T,e)]. In section |IIIFj we show 
that this critical point belongs to the Ising universality 
class. 

The shape of this phase diagram is identical to the one 
computed in na for a phenomenological hydrodynamic 
description of self-propelled particles with polar align¬ 
ment. The comparison with other microscopic models in 
the literature is however hard to make since there seems 
to be very few studies in the (e, po) plane, probably be¬ 
cause very few models admit a well-defined zero velocity 
limit. 


C. Nucleation vs spinodal decomposition 

As for an equilibrium liquid-gas transition, the coex¬ 
istence lines pg{T,e) and p£{T,e) are complemented by 
spinodal lines ^g{T^ e) and (p^(T, e) that mark the limit of 
linear stability of the homogeneous gas and liquid phases. 



FIG. 7. Successive snapshots following quenches from homo¬ 
geneous gas and liquid phases inside and outside the spinodal 
region. Parameters: D = 1,£ = 0.9,/3 = 1.8, system sizes 
400x400 and 1000x1000 for the quenches from the gas and 
liquid phases. From top to bottom, po = 1.84, 3, 3, 4.7. See 
Supplementary Movies in [49j . 


respectively. While pg and p£ are easily measured in sim¬ 
ulations, (fg and (p£ are much harder to access numerically 
at non-zero temperature: When the system is in the co¬ 
existence region but outside the putative spinodal lines, 
the homogeneous phases are metastable and finite fluctu¬ 
ations make the system phase-separate. The closer to the 
spinodal line, the faster this nucleation occurs and it is 
then difficult to pinpoint precisely the transition from a 
‘fast’ nucleation to a spinodal decomposition. Neverthe¬ 
less, the differences between the coexistence and spinodal 
regions are clearly seen when, starting from a homoge¬ 
neous phase, one quenches the system in the coexistence 
region but relatively far away from the spinodal lines. 

Quenching outside the spinodal region, the homoge¬ 
neous phases are metastable. The closer to the binodals, 
the longer it takes for a liquid (resp. gas) domain to be 
nucleated in the gas (resp. liquid) background. The con¬ 
vergence to the phase-separated steady-state then results 
from the coarsening of this domain. 

Quenching inside the spinodal region, the different 
symmetries between gas and liquid result in different 
spinodal decomposition dynamics when starting from or¬ 
dered and disordered phases. Starting from a disordered 
gas, the linear instability almost immediately results in 
the formation of an extensive number of small clusters of 
negative and positive spins. The coarsening then stems 
from the merging of these clusters, until a single, macro¬ 
scopic domain remains. The late stage of the coarsening 
is thus dominated by the long-lived competition between 
a small number of right- and left-moving macroscopic 
clusters. Their shapes (see Fig. are reminiscent of 
the counter-propagating arrays of bands reported in m, 
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FIG. 8. Snashots in the late stage of coarsening taken from 
the same simulation as the first row of fig. at time t = 
283000 (left) and t = 310000 (right). Parameters: D — — 

0.9,/3 = 1.8, yoo = 3, system sizes 400x400. 


where it was suggested, using deterministic simulations 
of the Boltzmann equation derived for kinetic flocking 
models, that such profiles could constitute a new phase of 
flocking models. In our simulations, we always observed a 
coarsening process leading to a single band, which seems 
to indicate that the apparent stability of these solutions 
in m could be due to the lack of fluctuation terms. It 
would nevertheless be interesting to make a more detailed 
study of the coarsening dynamics to see if these alternat¬ 
ing bands could indeed form a stable phase (for instance 
at low temperatures, where the coarsening seems to be¬ 
come slower and slower). 

Starting from the ordered phase, the linear instabil¬ 
ity results in many liquid domains which all move in the 
same direction. The coarsening then results from the 
collision of liquid bands that move in the same direc¬ 
tion, but with slightly different speeds. See Fig. 0and SI 
movies [l^ for examples of these four possible dynamics. 


D. Hysteresis loops 

Another similarity with a liquid-gas transition is the 
presence of hysteresis loops obtained by varying slowly 
the density at constant p and 5 in finite-size systems. 
Such loops are shown in the left panel of Fig.[^ where the 
liquid fraction ^ is reported as the density is continuously 
ramped up and down. To measure ^ numerically, a first 
strategy, followed in [22], is to compute average density 
profiles at fixed po? as shown in the left panel of Fig. [^ 
and use an arbitrary density treshold between pg and pi 
to associate each site to the gas or liquid regions. 

Since the interfaces between gas and liquids are not 
perfectly straight, this is slightly artefactual for finite- 
size systems. Here we decided to measure ^ numerically 
through: 


^num. = 




E 


rui 


( 4 ) 


where mi is the magnetization of the plateau in the liq¬ 
uid part of the profile, {mi is independent of po as long 
as the system is phase-separated and corresponds to the 


FIG. 9. Left: Evolution of the liquid fraction cj) upon 
changing continuously po. Large jumps in 4> correspond to 
the nucleation of bands in meta-stable homogeneous pro¬ 
files while small jumps are finite-size effects due to the finite 
width of the interfaces connecting the gas and liquid regions. 
The density is increased by ^po = 0.02 every At = 2000. 
Lx X Ly = 800 X 100, yd = 2, £ = 0.9, D = 1. Right Evolu¬ 
tion of the magnetisations per particles and per sites as the 
density is varied. The linear scaling typical of liquid-gas phase 
transition is seen using tul = Mjl?. 


magnetization of a uniform liquid phase at the coexis¬ 
tence density p^.) The results are very similar to those 
obtained in [22] but first 0 is faster to measure and 
second it does not rely on an arbitrary density treshold. 

Starting in the gas phase and increasing the density, 
the system remains disordered, with a liquid fraction 
0 = 0, until a band of liquid is nucleated, at which point 
0 jumps to a finite value. Increasing again po, the liquid 
region widens until the two interfaces between gas and 
liquid almost touch and the liquid phase almost fills the 
system. At that point, the system jumps to a homoge¬ 
neous liquid phase with 0=1. 

Upon decreasing the density, a similar scenario occurs: 
A homogeneous liquid becomes metastable as the coexis¬ 
tence line is crossed. As the density keeps decreasing, the 
system thus remains in a liquid state with ^ = 1 until 
a nucleation event brings it to a phase-separated profile. 
The liquid region then shrinks until its boundaries al¬ 
most touch and a second discontinuity of ^ occurs as the 
system jumps into a homogeneous gas phase. 


E. Order parameter and finite-size scaling 

The liquid-gas transition picture suggests different 
finite-size scaling and order parameter than those asso¬ 
ciated to magnetic phase transitions previously used to 
study flocking models. Most studies [mini EH indeed 
relied on the mean magnetization per particle 

= ( 5 ) 

i 

rather than the mean magnetization per unit area 

niL = A— W TOj = porriN (6) 

LxLy . 
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FIG. 10. Hysteresis loops for system sizes 200x100, 400x100 
and 800x100. For each system, the density is increased by 
(5po = 0.02 every At = 2000. T = 0.5, £ = 0.9. 

For models like the Vicsek model, the former is noth¬ 
ing but the polarisation iiitv = P while the latter is 
related to the total momentum mL = poV/v. In the 
phase-separated region, both can be related to the liquid 
fraction ^ through Eq. © 

r 1-Po/PO 

rriN = -^^L^Lyini = mi -- ( 7 ) 

A pi - pg 

mL = ^mi = mi——— ( 8 ) 

Pi- Pg 

The simple linear scaling of mL with pQ — pg is replaced 
by a non-linear dependence of rriAr with po, as shown in 
Fig. in right panel. An apparently inoccent change of 
the normalization used to make the magnetisation M = 
mi intensive can thus turn the simple affine scaling 
of with po, typical of a liquid-gas transition, into 
the non-linear dependence of ttt-at that could make one 
mistake the transition for a critical one. 

Let us now go back to the hysteresis loops and discuss 
their finite size scaling. As shown in figure [T q| the discon¬ 
tinuities of the liquid fraction get closer and closer to the 
binodals pg and pi as the system size increases, leading to 
vanishingly small hysteresis loops in the thermodynamic 
limit. 

Consider first the transition from gas to phase- 
separated profiles. The liquid fraction exhibits two dif¬ 
ferent discontinuities when the density is decreased or 
increased, due to two different effects. As the density 
is decreased, phase-separated profiles cannot be main¬ 
tained arbitrarily close to pg. There is indeed a critical 
nucleus, which roughly amounts to two connected domain 
walls, as can be seen in Fig. for po = 1-2. As shown 
on Fig. (left), this critical nucleus Lc is independant 
of the system size. If the excess mass LxLy{pQ — pg) 
is smaller than a critical value (pcLy^ this critical nu¬ 
cleus cannot be accomodated in the system, which thus 
falls into the gas phase. As the system size increases, 
the minimal density to observe phase-separated profiles 
Po = Pg -f- Pc/{Lx) thus converges to pg as Lx increases 
and phase-separated profiles are seen closer and closer to 
the binodal. The second discontinuity, met upon increas¬ 
ing the density, corresponds to the nucleation of a liquid 
band of width in a gaseous background. Since can 



FIG. 11. Left: Divergence of the critical nucleus Lc when 
approaching the critical points /S ^ 1 and £ ^ 0. To measure 
Lc, we started in the phase-separated state and decreased con¬ 
tinuously the density (the errorbars correspond to the density 
step used) to record the density pm at which the liquid band 
disappears. Lc is then defined by Lpm = Lpg + Lc{pi — pg), 
as the length of a band at density pi that can be made with 
the excess density pm — Pg- Right: Variation of the critical 
nucleus iwht L showing that, within numerical errors, it does 
not depend on system size. Parameters: D = 1, e = 0.9, 
13 = 1.9 


be anything between Lc and Lx, increasing the system 
size at fixed density should decrease the mean time until 
nucleation of such bands, thanks to an entropic contri¬ 
bution due to the number of places where the bands can 
be nucleated. As shown in this is indeed the case 
and the transition to phase-separated profiles thus also 
happens closer and closer to the binodal pg. 

The same line of reasoning can be used to understand 
the scaling of the second hysteresis window, close to pi. 
Thus, in the thermodynamic limit, all discontinuities dis¬ 
appear and the liquid fraction varies continuously from 
(j) = f) dX pQ = pg to (j) = 1 dX pQ = pi, as for an equi¬ 
librium liquid-gas transition in the canonical ensemble. 
Note that the width of the critical nucleus diverges as 
one gets closer and closer to the critical points {e = 0 
or (3 = 1), as shown in the right panel of Fig. Ku This 
could explain why some studies of the Vicsek model in 
the small velocity region claim to find a critical transi¬ 
tion [33]: as one gets closer and closer to the zero speed 
limit, the system-size above which one can correctly ob¬ 
serve the discontinuous nature of the transition diverges. 


F. The e = 0 critical point 

While the (3 = 1, pc = oo, critical point is out of 
reach numerically, the study of the 5 = 0 critical point 
is accessible. At 5 = 0, there is no self-propulsion and 
the phase transition is of a completely different nature 
from the liquid-gas transition described above. As we 
show below, despite the dynamics being non-equilibrium, 
it turns out to be a standard critical phase transition 
belonging to the Ising universality class. 

We studied this critical point using a finite-size scaling 
standard for magnetic systems at criticality [3l] . We thus 
consider the magnetization m^ G [0,1]. In equilibrium. 


















FIG. 12. Left: Binder cumulant G{p) from which we find the critical density p* — 2.798 zb 0.002. Other three figures: data 
collapse on the universal scaling functions iGi, and Fg (defined in the text) when the data is rescaled with the 2d Ising 

exponent values jS — 1/8, 7 = 7/4 and ly — 1. t — p^')!p^. Parameters: D — e — 0.9, (3 = 1.9. 


around ferromagnetic critical points, the order parame¬ 
ter, susceptibility and Binder cumulant G = 1 — 
are known [3l] to obey the finite-size scaling relations 


(m) = (9) 

X = L\{m?) - {mf) = (10) 

G = FcitL^^n (11) 

where t = L^^^{p — p*)/p* is the rescaled distance to the 
critical density p*. F^ and Fq are universal scaling 
functions and /3, 7 and v the usual critical exponents. 

We used the fact that G{t = 0) is independent of L to 
find the critical density, which is thus the density where 
all the curves G{t) for different system sizes intersect 
(Fig. left). We found that the value at the cross¬ 
ing point is the same universal value G{t = 0) 0.61 as 

in the 2d Ising model [35]. A very neat data collapse is 
further observed for the critical exponents of the 2d Ising 
model [3 = 1/8, 7 = 7/4 and u = 1 (see Fig. [^. We 
thus conclude that the critical point at 5 = 0 is indeed in 
the Ising universality class. 

Note that a direct evalutation of the critical exponents 
is much harder than for the equilibrium Ising model. 
Here, the dynamics is fixed so one cannot use alterna¬ 
tive dynamics like cluster algorithms to circumvent the 
problem of critical slowing down. 


G. Number fluctuations 

In most flocking models the homogeneous ordered 
phase exhibits giant density fluctuations [13 [13 EHl 133 
|37|. These are quantified by measuring number fluctu¬ 
ations, i.e by counting the number of particles n{i) in 
boxes of increasing sizes £ < L and computing its root 
mean square An{£). When the correlation length jC is 
finite, a box of size £^ C can be divided in (^/G)^ inde¬ 
pendant boxes. The total number of particles in the large 
box is then the sum of independent identically distributed 
random variables; The central limit theorem applies and 


the probability distribution of n{£) tends to a Gaussian. 
This yields the “normal” scaling An ^ On the 

contrary, one finds in the Vicsek model the anomalous 
scaling An ^ n^-^ [T3] . 

In the Active Ising model the number fluctuations are 
found to be normal in the liquid and gas phases, where 
An ^ and trivially ‘giant’ in the phase-separated 

regime where An ^ n (see Fig. 13). 

Note that the scaling An ^ n is a simple consequence 
of phase-separation and one should thus distinguish this 
scaling from the ‘anomalous’ scaling of the Vicsek model, 
which is a signature of long-range correlations. Let us 
consider a system at liquid fraction 0 that is large enough 
that we can find a range of box sizes £ such that: 1) £ <^ L 
so that we can neglect the contribution of the interfaces 
(a box is either in the liquid or the gas phase); 2) £ is 
large enough that n{£) takes only two possible values and 
we can neglect the fluctuations around these two values. 
With these assumptions. 


P(n) ~ ^(5(n - pef) + (1 - (())<5(n - pgi'^) (12) 

where pi and pg are the densities in the gas and liquid 
domains. Then one finds 

(n) = {<t)pt + (1 - (fypg)^^ = (13) 

An = V(n2) - (n)2 = 7<(,(1 - (14) 

Po 

which is a simple hand-waving explanation of the scal¬ 
ing observed in the coexistence region of the active 
Ising model, as well as in other phase-separating sys¬ 
tems 


IV. HYDRODYNAMIC DESCRIPTION OF THE 
ACTIVE ISING MODEL 

In this section we derive and analyze a continuous de¬ 
scription of the AIM based on two coupled partial differ¬ 
ential equations accounting for the spatio-temporal evo¬ 
lutions of the density and magnetization fields. 
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FIG. 13. Number fluctuations in the three different phases: 
gas (red), liquid (green) and at coexistence (blue, upper line). 
n is the number of particles in boxes of size i and An its root 
mean square. D = 1, 5 = 0.9, L = 400, po = 5 


We first show in section IIV Al that a standard mean- 
field treatment wrongly predicts a continuous transition 
between the disordered gas and the ordered liquid. In 
section |IVB| we show that local fluctuations, which are 
neglected in the mean-field approximations, are neces¬ 
sary to correctly account for the physics of the system 
when the density is finite and 5 7 ^ 0. We show in particu¬ 
lar that as soon as the density is finite, fluctuations make 
the transition first order. We then use our hydrodynamic 
description in section |V| to study the inhomogeneous pro¬ 
files. 


A. Mean-field equations 

The simplest way to account analytically for a non¬ 
equilibrium lattice gas is probably to derive mean-field 
equations. These are known to be quantitatively wrong, 
but they often capture phase diagrams correctly gQlSI]. 

Their derivations follow a standard procedure which 
can be applied to the AIM and which, for simplicity, we 
first present in ID. Starting from the master equation, 
one first derives the time-evolution of the mean number 
of ±1 spins on site i 

{hf) = D{1 ± e){nU) + D{1 t e){n%^) - 2D{nf) 

± (n“ exp(/3—)) =F (N exp(-/3—)) (15) 

Pi Pi 

which can then be rewritten for the density and magneti¬ 
sation 


{pi) = D{{pi+i) + {pi-i) - 2{pi)) - De{{rrn+x) - {rrn-x)) 

(16) 

(mi) = D{{mi+i) + (mi_i) - 2{mi)) - De{{pi+i) - (pi-i)) 

771 ■ 771 ■ 

+ 2{pi sinh(^^)) - 2(mi cosh(^^)) (17) 

Pi Pi 

One can then take a continuum limit using the rescaled 
variable x = ijL G [0,1], D = Djl?, v = 2DelL and 
use the Taylor expansion pi±i = p{x) ± L~^dxp{x) + 
L~‘^dxxp{x)/2. We then obtain equations for the contin¬ 
uum fields p(x), m(x), which are assumed to smoothly 


interpolate the discrete occupancies p^, 


dt{p) = Ddxx{p) - vdx{m) 


(18) 


dt{m) = Ddxxi'^) ~ ^dx{p) + ^ 2 psinh — 2 mcosh 

^ m 


In higher dimension, the sole difference is that the dif¬ 
fusive terms become DA{p) and DA{m) whereas the v 
terms still involve solely dx since the hopping is biased 
only horizontally. 

In practice, to compare microscopic simulations and 
hydrodynamic theories it is often easier not to rescale 
space and use a continuous variable x = Lx e [ 0 , 1 /] (and 
hence D = DL^ and v = Lv = 2De). Macroscopic and 
microscopic transport parameters are then expressed in 
the same units and equations and are then valid, 
without the tilde variables. This is what we use in the 
following. 

Equations and (19) are exact; they couple the 
first moments^) and to higher moments through 
the averages of the hyperbolic sine and cosine functions. 
Following the standard procedure established for equi¬ 
librium ferromagnetic models, we then make two ap¬ 
proximations. First, we take a mean-field approximation 
by replacing (/(p, m)) by /((p), (m)), for any function 
/. (We then drop the (...) notation for clarity.) This 
amounts to neglecting both the correlations between den¬ 
sity and magnetisation and their fluctuations. Second, 
we expand the hyperbolic functions in power series, up 
to w? jp^. This further restricts our description to the 
case where m p. We then arrive at the mean-field 
equations 


p = DAp — vdxTn 


( 20 ) 




rfi = DAm — vdxp + 2m(/3 — 1 ) — ( 21 ) 

p^ 

where a = (3^{1 — P/3). (For /3 > 3, one should expand 
to higher order to obtain a stabilizing term.) 

Let us consider the various terms appearing in the 
mean-field equations. The first terms on the r.h.s of (20) 


and ( 21 ) are diffusion terms arising from the stochastic 
particle hopping. Let us stress that these terms do not 
depend on the bias e and are thus present even in the 
totally asymmetric case 5 = 1; they do not rely on the 
possibility for +1 and —1 particles to hop backwards and 
forwards, respectively. The second terms, proportional 
to are due to the bias. Their physical origin is ex¬ 
plained in Fig. where we show how positive gradients 
in m or p yield negative contributions to p or m, respec¬ 
tively. Finally, the last two terms in ( 21 ) stem from the 
ferromagnetic interaction and, apart from the p^ depen¬ 
dence of the last term, are typical of Landau mean- 
field theory. Note that the alignment terms are the only 
non-linear ones and thus the only terms for which the 
mean-field approximation is actually an approximation. 
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FIG. 14. Schematic account for the role of density and mag¬ 
netisation gradients in the mean-fie Id equations. Left: Ini¬ 
tially, m = 0 and Vp > 0 around site i. Once plus particles 
jump to the right and minus particles to the left, the density 
in site i is unchanged but rrn has decreased. Right: Ini¬ 
tially, p is constant and Vm < 0 around site i. Once particles 
have jumped, the magnetisation of site i is unchanged but the 
density has decreased. 

The mean-field equations always accept the trivial ho¬ 
mogeneous solution 

p{x)=po, m{x) = 0, (22) 


FIG. 15. Linear stability of homogeneous profiles in the naive 
(left) and refined (right) mean-field models. Plain (resp. 
dashed) lines denote stable (resp. unstable) solutions. In 
the RMFM, for /3 < 1, only the homogeneous profile exists 
and is stable at all densities. 

When computing, for instance, the first non-linear 
term neglecting correlations between m and p 

leads to 

(toVP) = (24) 


which is linearly stable for /3 < 1. As soon as /3 > 1, two 
ordered homogeneous solutions appear, 

P = P0, m = ±po\ — —(23) 

which are linearly stable (see the left panel of Fig. [T^. 
Therefore, at the mean-field level, a linearly stable ho¬ 
mogeneous solution exists for all (/3, po). Furthermore, 
integrating numerically Eqs. ([2Q| ) and ( [21] ) starting from 
different initial conditions |47j/ the system always re¬ 
laxes to a homogeneous solution and inhomogeneous pro¬ 
files are never observed. Hence, the mean-field equa¬ 
tions predict a continuous transition between homoge¬ 
neous disordered and ordered profiles at /3 = 1, just as 
for the Weiss ferromagnet [l2]. The phase diagram is 
simply split between a high-temperature disordered ho¬ 
mogeneous phase, for T > 1, and a low temperature or¬ 
dered homogeneous phase, for T < 1. This mean-field 
approach thus completely misses the phenomenology of 
the microscopic model; It cannot explain the existence of 
phase-separated profiles and yield a phase diagram corre¬ 
sponding to a single (continuous) transition line at T = 1, 
in contradiction to the crescent shape observed in the mi¬ 
croscopic model (see Fig. [^. 


B. Going beyond the mean-field approximation 


Previous coarse-graining approaches of flocking mod¬ 
els [m na 113 E] often relied on neglecting correlations 
by factorizing probability distributions. For example, in 
the Boltzmann-Ginzburg-Landau approach of Bertin et 
al. ca the two-particle probability distribution is re¬ 
placed by the product of one-particle distributions. In 
our case, to derive the mean-field equations (20) and (21) 
we made an even cruder approximation. 


We went one step further, completely discarding fluctu¬ 
ations and replaced (1/p^), (m^) by l/(p)^, (m)^. As 
we show below, these fluctuations are crucial to account 
qualitatively for the physics of the AIM. 

The dynamical equations ( [l8| and (B^ on the first 
moments predict how (p(x,t)yand (m^t)) evolve in 
time, given an initial distribution 

P[/9, m] = 5{p{x) - po)5{m{x) - toq). (25) 


The mean-field approximation then amounts to compute 
the averages of hyperbolic functions in (19) by assum¬ 


ing that, as time goes on, V remains a product of Dirac 
functions: 


vip, m;x,t\po, mo] = 5{p{x, t) - p{x, t)) 

X 6{m{x, t) — ffi{x^ t)) 


(26) 


where p(x, t) and ffi{x^) are solutions of the mean-field 
equations ( [^ and ( j^ . In practice this means that re¬ 
peatedly simulating the microscopic model starting from 
an initial distribution ( |2^ always yields the exact same 
values p{x,t) = p{x^t) and m{x,t) = m{x,t). A bet¬ 
ter description should allow both p and m to fluctuate 
around their mean values as well as account for their cor¬ 
relations. 

We can thus improve our approximation by replac¬ 
ing the dirac functions in ( |26| ) by Gaussians of variance 
cr^(x,t) and cr^(x,t). This still neglects correlations be¬ 
tween p and m but allows for (small) fluctuations around 
their mean. Note that the only approximation made in 
the derivation of the mean-field equations occured at the 
level of the alignment terms. Since each site of the AIM 
is a fully connected Ising model, it is reasonnable to as¬ 
sume that in the large density limit, mean-field should be 
correct. We thus assume that our corrections to mean- 
field should be small in the high density regions, where 
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it is also reasonnable to assume that the variance of 
p(x, t) and m(x, t) are proportional to p\ = app and 
= ^mp where ap and am are functions of /3 and v 
only. 

The probability to observe given values of p(x, t) and 
m(x, t) is then assumed to be 


V[p, m;x,t\ po, ^o] = A/'(p - p, app)M{m - m, amp) 

(27) 

where /\/27rcr^ is the normal distri¬ 

bution. 


Under these assumptions, the alignment term in (19) 


can still be computed analytically; We show in ap¬ 
pendix]^ that, at leading order in a m/p expansion, the 
correction to mean-field reads 

pm \ 


/ 2 p sinh — 2m cosh ^- / 

\ P P / 


2ip — 1 — -)m — a^TT 

P P^ 

(28) 

where r = Zaaml‘2^ is a positive function of p. Intu¬ 
itively, the fluctuations “renormalize” the transition tem¬ 
perature 


Pt{p) = 1 +!: = ^ 

p p 


(29) 


In principle, one could expand pt to higher order to ob¬ 
tain a better and better approximation. The correc¬ 
tion (29) however suffices to account qualitatively for the 


most salient features of the microscopic model and we 
will thus stop our expansion at this order. Furthermore, 
extending (28) to higher orders does not suffice to provide 


quantitative agreement between microscopic simulations 
of the AIM and the “corrected” mean-field equations, 
probably because we still neglect correlations between 
p and m. More details are provided in appendix for 
the interested reader. 


C. Refined Mean-Field Model 

The correction to mean-field derived in the previous 
section can thus be seen as a finite-density correction 
to the transition temperature Pt, which converges to its 
mean-field value p^^ = 1 as p ^ oo. As was already 
recognized in previous studies [13 na 123113], the density- 
dependence of Pt is the key ingredient to describe phase 
separation at the level of hydrodynamic equations. With 
this correction, we obtain a refined mean-field model 
(RMFM) 

p = D/S.p — vdxm (30) 

Q 

T m 

m = DAm — vdxP + 2ip — 1-)m — a^r (31) 

p p^ 

which we now study. 

The linear stability analysis of homogeneous solutions 
strongly differs from the mean-field case. For P > 1 the 
disordered profile is stable for po G [0,(p^(/3)] where 

</^s(/3) =^/(/5-1) (32) 



FIG. 16. Phase diagrams in the RMFM. The lines pg and 
Pt are the spinodal lines denoting the limit of linear stability 
of homogeneous profiles. The lines pg and pt are coexistence 
lines that limit the domain of existence of phase-separated 
profiles. Top row: temperature/density ensemble. The right 
plot is a zoom of the region delimited by the grey rectangle. 
D — r — V — 1. Bottom left: velocity/density ensemble. 
D — r — P — l.b. Bottom right: 2d snapshot of the 
density field in the phase coexistence region. Its position in 
the phase diagrams is indicated by the grey squares. D — r — 
V = 1, P — l.b and Po = 2.1. 


The homogeneous ordered solutions 


p{x) = Po, 


m{x) = mo = ±po 



- 2 - 


poa 


(33) 


exist for all po > pg, but are only stable for po > > 

Pg (see Fig. [^. The explicit expression of pg can be 
found using a standard linear stability analysis, detailed 
in Appendix [B] 




v^/a + 8D{P — 1 )^) + v‘^K, -h 8Da{P — 1) 


2v^K + SDa{P -1) 


where k = 2 a — 2p. 

Close to the critical point at P 


(34) 


1 , 


+ X-h 0{p — 1) (35) 


so that Pt and Pg both diverge, while their difference 
remains constant. Close to the u = 0 critical point, we 
obtain 




rv 


VSDa{p - 1 ) 


+ 0(u') 


(36) 


so that Pt Pg when u ^ 0. 

The homogeneous solutions are linearly unstable in the 
density range [pg,pt]. Simulating the RMFM [37] for 
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FIG. 17. Hysteresis loops in the RMFM. Left: Density pro¬ 
files along the loop as po is varied. Right: Evolution of the 
liquid fraction 0 upon changing continuously po- Parameters: 
^ = 1.5, D = i; = r = l, L = 800. 


such densities yield phase-separated profiles similar to 
those seen in the AIM, with macroscopic liquid bands 
travelling in a disordered gas background (see bottom- 
right panel of fig. 16). The densities in the gas and liquid 
parts of the profiles remain constant as po is varied; they 
thus give access to the coexistence lines pg and pi. 

The phase diagrams of the RMFM in the tempera¬ 
ture/density and velocity/density ensembles shown in 
Fig. are qualitatively similar to those of the AIM, 
with an asymptote at T = 1 when po —> oc in the (T, po) 
plane, and a critical point at u = 0 in the (u, po) plane. 
As before, the coexistence lines pg and pi delimit the 
domain of existence of phase-separated solutions; they 
can now be complemented by the spinodals cpg and cpi 
which mark the loss of linear stability of homogeneous 
disordered and order phases, respectively. 

The hysteresis loops observed in the RMFM (see 
Fig. are similar to those found in the microscopic 
modd^see Fig. and |^. Starting at low density in the 
gaseous phase and increasing density the system stays 
in the gas phase until it becomes unstable at po = 
where it phase separates. Increasing again the density, 
the liquid fraction increases linearly until the liquid al¬ 
most fills the system. As in the AIM, the finite widths of 
the interfaces set a minimum and a maximum size for a 
domain, hence preventing liquid bands from completely 
filling the system. This results in a discontinuous jump 
of the liquid fraction close to the binodals, whose height 
vanishes as the system size diverges (see Fig. right 
panel). The main difference with the hysteresis loops ob¬ 
served for the AIM is that, given the absence of noise in 
the RMFM, there is no nucleation and the system phase- 
separates only when the spinodal densities are reached. 


D. Control parameters 


To determine how many independent control parame¬ 
ters are needed to describe the behavior of the RMFM, 
we recast Eqs. (30) and (31) in dimensionless form. To 
do so, we first have to introduce back the rate 7 which 
appeared in the definition of the flipping rates § and 
that we have taken equal to one until now. Introducing 


the dimensionless variables and constants 


t = t/ 7 , X = A —x, p = rp, m = rrh^ 


(37) 


the refined mean-field equations become 

p = Ap — (38) 

• ^ 1 TTL^ 

m = Am — vdxp + [2(/3 — 1 -)m — (39) 

^ p p^ 


Since a is a function of /3, there are only two external 
dimensionless control parameters: u is a Peclet number 
comparing the advection speed v and the diffusivity D at 
the length scale u /7 travelled by a particle between two 
spin flips; /3 which controls the ordering of the system [ 38 ] 
The average density, which sets an external constraint on 
the system, constitutes a third independent parameter. 
Our phase diagrams shown in Fig. thus sample all the 
relevant parameters of the RMFM. 


V. INHOMOGENEOUS BAND PROFILES 

In the previous section we have shown how one can 
build a refined mean-field model by taking into account 
the local fluctuations of magnetisations and densities. 
Numerical simulations of the RMFM exhibit a phe¬ 
nomenology akin to that of the microscopic AIM, con¬ 
firming the liquid-gas picture of the phase transition. We 
now focus on the inhomogeneous profiles and show ana¬ 
lytically that the RMFM accounts for their shapes and 
speeds when (3^1. Furthermore, the RMFM also cor¬ 
rectly predicts the scaling of the width of the critical 
bands in the vicinity of the critical points (3 = 1 and 

V = 0. 


A. Propagative solutions 


Let us reduce Eqs. (30) and (31) to a single ordinary 


differential equation. To do so, we first introduce a new 
coordinate z = x — ct comoving with the liquid band 
at an unknown speed c. In this comoving frame, the 
stationnary solutions of the RMFM satisfy 


Dp" + cp' — vm' = 0 
Dm" + cm' — vp' + 2{f3 — 1 


(40) 

Q 

r TD 

-)m-aAr=0 (41) 

p p^ 


The RMFM is a finite-density correction to the p = 
00 mean-field limit and should thus work best for large 
densities. As we can see on the phase diagram shown 
in Fig the densities pi and pg diverge as /3 ^ 1, as 
do ifg and pi (see Eqs. (34) and (35)). Furthermore, 
one can check that pi — pg remains finite in this limit, 
as does mi (see Fig. 18). Close to /3 = 1 , we can thus 
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expand Eq. (41) in power of e = mj^pg ~ dpjpg^ where 
6p={p- <Pg)7to get 


0 = Dm" + cm' — vSp' + 


2rmSp m 


‘fl 


- a— (42) 


Besides, Eq. (40) can be solved iteratively to obtain 


p{z) in terms of rn\z) and its derivatives 


p[z) = Pg + -m[z) + - 
c c 


oo 

i:( 


k=l 


D\^ d^m{z) 
c ) dz^ 


(43) 


where pg is an integration constant that equals the den¬ 
sity in the gas phase at coexistence, since p{z) = pg where 
m = 0. Again, the RMEM should work best close to the 
critical points, where the width of band fronts diverge 
(see Eig. 11), we can thus expect the development (43) 
to rapidly converge in this limit and retain only 


/ N N /// N 

p[z) = Pg ^ -m{z) - 


(44) 


At second order in e, Eqs (42) and (44) then reduces 


to 


Dm" + (ao — aim)m' — bim + h 2 m^ — = 0 (45) 


where we have introduced the positive constants 



FIG. 18. Left: The magnetization rrii in the liquid band and 
Pi — Pg at phase-coexistence, measured in the microscopic 
simulations, converge to the same constant when ^ ^ 1 as 
predicted by the analytical solution. Parameters: D — e — 
0.9, L — 400 for the microscopic simulations, r — v — D — 

L = 400 for the RMFM. Right: velocity c of a liquid band 
propagating in a gas background. As /3 ^ 1, c u i n the 
microscopic model, in Id simulation of the RMFM (31), and 
in the analytical solution. 


but straightforward algebra then gives 
Sr^D 


Arv 

3(ac 


c = v(^l 

i 

mi 

Pg ~ ^9 ~ 

^a/d 


3av‘^p‘l y 


4rv^ 


(48) 


9ac^ 

C7 




1±1 1 + 


47 + 

3 


'y:^av^(Pg 

2Dr2 


7 ; 7 ; 4 DllV 


bi = 2r 


Pg Pg 

pI 


bo = 


2rv 

9 ’ 


u “ 

h = ^ 


(46) 

We then look for propagating solutions made of two 
fronts, connecting an ordered liquid band at mi to a 
disordered gas background at pg^ mg = 0. Precisely, we 
look for propagating fronts given by: 


m{z) = ^ [1 + tanh(A: 2 :)] 


(47) 


To describe phase-separated domains, we need two front 
solutions, an ascending front ma{z) with ka > 0 and 
a descending front md{z) with kd < 0, with the same 
speed c, density pg and magnetization m^. Since the 
term (uq — aim)m' breaks the symmetry of the equations 
under (m^c) {—m^—c) the fore and rear fronts need 

not be the same, so that \ka\ 7 ^ \kd\ in general. 

The complete solution, specified by (c, pg^ mi, Pa/d) 
can be obtained by injecting the Ansatz (ITl) into 
Eq. (45). Using the equality tanh'(A:x) = /c —/c tanh^(/cx), 
the l.h.s. of Eq. (45) then yields a third order polynomial 


in tanh(/cx) whose coefficients all have to vanish. Tedious 


where 7 ± = 1 ± ^. 

The solution is thus completely determined, the den¬ 
sity and magnetization profiles being given by Eq. (47) 
and (44) respectively. 


B. Close to the /3 = 1 critical point 

At leading orders when /3 ^ 1 , the propagating fronts 
are characterized by 


Pg ~ Pg 


4r 


P£ — Pg P 




(49) 


me = 

3a 


ka = 


4r 

3a 

p' 


3V Da 


6va 


c = V -\- 
kd= - 


2D{p-l)‘^ ^ 

3va 

P-I (^- 1)2 


3V Da 


6va 


Some comments are in order. Eirst, the two coexis¬ 
tence lines pg and pi diverge as /3 ^ 1 , as do the spinodals 
Pg and pi^ while their difference and the magnetization 
mi converge to finite constants. This behavior, which 
is in line with simulations of the microscopic model (see 
Eig. p! 8 |), legitimates the expansion of in powers of 
mjpg and dpjpg. 
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Then, we can check the validity of the gradient ex¬ 
pansion by comparing two successive terms in Eq. (43). 
When /3 ^ 1, we have 


K c ) 

/ D \ ^ d^m 
\ c ) dy^ 


^^a/d 


(^- 1 ) 


(50) 


SO that our approximation becomes exact when /3 ^ 1. 

The front solutions account for a number of interesting 
features of the propagating liquid bands. First, the front 
speed c is generally larger than the maximal mean 
speed of a single spin. This may seem surprising un¬ 
til one realizes that the front propagation is due both 
to the spins in the liquid band hopping forward and to 
the “conversion” of disordered sites into ordered ones at 
the level of the fore front. There is thus a FKPP-like 
contribution m to the speed of a band, which allows c 
to be larger than v. Interestingly, despite the approx¬ 
imations made in deriving the RMFM, the behavior of 
c/v diS /3 ^ 1 coincides exactly with what is observed in 
microscopic simulations of the AIM (see Fig. 18). 

Regarding the propagating fronts, the analytical solu¬ 
tion predicts \ka\ < \kd\, i.e., that the descending (fore) 
front is steeper than the ascending (rear) front. The 
asymmetric term being subleading as /3 ^ 1 , the fore and 
rear fronts become more and more symmetric as /3 ^ 1 . 
This is consistent with the microscopic model: In Fig.p!^ 
we show that the fronts are well described by two sym¬ 
metric tanh functions close to /3 = 1. As the temperature 
decreases, the fronts first remain well approximated by 
hyperbolic tangents, but with different widths ka 7 ^ kd^ 
before their functional form deviates from the tanh solu¬ 


tion (see Fig. 19). 

Let us now be slighlty more quantitative and compare 
the scalings of the front widths in the AIM with the pre¬ 
diction of our analytical solution (49). In the microscopic 


model, we fitted the fronts of phase-separated profiles 
by the hyperbolic tangent solutions (47) to extract their 


width. Although data is hard to obtain close to critical 
points, because m/p 0 , the measures are consistent 
with the analytical predictions. As shown in Fig. 

^a/d ^ (/^ “ 1) when P ^ 1. One can also see that in this 
limit the two fronts become symmetric, i.e., ka kd- 
The size of the interfaces, inversely proportional to ka/d^ 
can be linked to the size of the critical nucleus. As ex¬ 
plained in sec. IIID[ a liquid domain can form only if 
the excess number of particles with respect to the gas is 
sufficient to create a band of minimal size Lc. As a first 
approximation, this minimal size is set by the size of the 
interfaces so that we expect Lc ^ 1/ka ^ l/kd^ Indeed, 
the same scalings are observed for Lc as for l/kajd as 
shown in Fig. [TO 


C. Close to the = 0 critical points 

While our approach was derived to work close to the 
critical point at /3 = I, the front solution still predicts 


many correct scalings close to the v = 0 critical points. 
There, the propagating bands are characterized by 


Pg ~ Pg 

pe = + 


a/2 ) 


rv 


3(13 - l)V^Da 
\/8rv 

3(/3 - l)V3Da 

A \ 1/4 


(51) 


"•“[ 27W iraP ^ 

V 12DV6^ 

V udVq^ 


C = 


kn, - 


kd = 


Again, the two coexisting densities merge with the 
spinodal lines at v = 0 while the magnetization in the liq¬ 
uid vanishes, hence justifying the expansion of Eq. 
in powers oimjpg and 6p/ipg. While gradients are again 
expected to vanish as 1 ; ^ 0 , the expansion of p in deriva¬ 
tives of m includes a diverging prefactor (D/c)^ ^ l/v^^‘^ 
at the k^^ order. The comparison of two successive terms 
in the expansion (44) then yields 


( 


D\k+1 d^+^r 






/ D \ ^ d^m 
\ c ) dy^ 


Pk 


a/d ' 


0 ( 1 ) 


(52) 


Thus, in this limit, the series may still converge but the 
ratios between two consecutive terms do not vanish as 
'T’ ^ 0 and we cannot completely neglect higher order 
gradients. Nevertheless, as shown in Fig. the an¬ 
alytical solution correctly predicts that the asymmetry 
between the fore and rear fronts does not disappear in 
the 'T’ ^ 0 limit. It also correctly predicts the scaling of 
the front widths ka/d ^ and thus the scaling of the 
critical nucleus in this regime. 

Beyond accounting for the shape of the phase diagram 
and the liquid-gas nature of the transition, the RMFM 
can thus correctly predict the shape of the band, their 
speed and the scaling of the critical nucleus in the vicin¬ 
ity of the critical points. In order to get a more quantita¬ 
tive agreement between the RMFM and the microscopic 
model, beyond the estimation of the unknown parameter 
r, one would probably needs to account for the corre¬ 
lations between m and p. Apart from quantitative cor¬ 
rections, these correlations however do not seem to play 
any role in controling the structure of the phase transi¬ 
tion and most features of the propagating bands. Inter¬ 
estingly, symmetric hyperbolic tangent front were also 
observed in hydrodynamic equations for self-propelled 
rods [29], even though in that case the domains are not 
moving. 
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FIG. 19. The fore and rear fronts of propagating bands become more asymmetric as T decreases. The shape of the fronts in the 
microscopic model (red curves) also deviate more and more from the analytical tanh solution (valid in the limit (5^1). Black 
dashed curves are fits of the rescaled fronts by expression ( [47| , where k is used as a htting parameter. Parameters: D — 
e — 0.9. Fronts are averaged over time and along the vertical direction. 



FIG. 20. Scaling of the front widths close to the critical points 
{3^1 (left) and i’ ^ 0 (right). The data is consistent with the 
predictions from the RMFM in these limits Eq. ( [4^ and (51) 
both for the scaling of ka/d and for the asymmetry between 
the ascending and descending fronts. £ = 0.9 (left), {3 — 1.^ 
(right) and D = 1. 


VI. ROBUSTNESS OF THE RESULTS 



0 200 400 600 800 1000 



FIG. 21. Active Ising model with closed boundary condi¬ 
tions. Top: snapshot of the density field. Bottom: space-time 
graph (averaged in the y-direction) showing the liquid domain 
bouncing back and forth in the box. Parameters: {3 — 1.9, 
po = 3, U = 1, £ = 0.9. See supplementary movies in 09]. 


Let us now discuss how the results presented in the pre¬ 
vious sections extend beyond our lattice-gas model with 
periodic boundary conditions. To do so, we consider the 
case of closed boundary conditions in section |VI A and 
study an off-lattice version of the AIM in section | VIB 

A. Closed boundary conditions 


phase-separation between a liquid domain and a gaseous 
disordered background (see Fig. top). When the liq¬ 
uid domain reaches a boundary, it accumulates close to 
the wall until its magnetisation flips, and crosses back the 
system in the other direction. This leads to the bouncing 
wave shown on Fig. (bottom), which is reminiscent of 
what is observed experimentally for the collective motion 
of colloidal rollers (see supplementary movies of 0 )- 


Since the ordered liquid domains always span the whole 
system in the vertical direction and propagates period¬ 
ically in the horizontal direction, one could think that 
their existence and stability relies on the use of periodic 
boundary conditions. To check this, we simulated the 
AIM in closed boxes. We tried different conditions at 
the edges of the box: When particles hit a wall, their 
spins were either flipped, randomized, or left unaltered. 

The same behavior was observed in all cases. First, one 
notice a small accumulation of particles close to the wall, 
which is typical of self-propelled particles 06]. Then, 
the system shows the same type of travelling bands as 
with periodic boundary conditions, with a macroscopic 


B. OfF-lattice version 

To show that the phenomenology of the AIM does not 
rely on the spatial discreteness of this lattice gas, we 
devised an off-lattice version of our model. To do so, we 
consider N particles in a continuous space of size LxXLy. 
Each particle carries a spin ±1, which flips at rate 

7T7 • 

W{s ^ —s) = exp(—/3 —-) (53) 

Pi 

where the local density pi and magnetization rrii are com¬ 
puted in disks of radius 1 . 
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FIG. 22. Phase diagram and phase-separated prohles for the 
off-lattice model showing the same behavior as the lattice 
model. Parameters: D — v — 1 and = 1.6 for the 
prohles. 


The position of the particle evolves according to the 
Langevin equation 


ti = Sive^^ + V2Dr] (54) 

where ri and Si are the position and spin of particle i and 
ry is a Gaussian white noise of unit variance. 

The phenomenology of this model is very similar to 
that of the AIM; Its phase diagram in the temperature- 
density ensemble shows the same three regions, with an 
asymptote at /3 = 1 as p ^ oo (Fig. left). As 
in the lattice model, only the liquid fraction changes 
when po is increased at hxed temperature as shown in 
Fig. (right). 


VII. DISCUSSION AND OUTLOOK 

In this paper we have characterized in detail the tran¬ 
sition to collective motion in the 2d active Ising model. 
For any temperature T < 1 and self-propulsion velocity 
'T’ > 0, there is a finite density range for which the sys¬ 
tem phase-separates into a polar liquid and a disordered 
gas. The densities at coexistence do not depend on T 
or V so that changing the average density only changes 
the liquid fraction. This is one of the many character¬ 
istics shared by the flocking transition of the AIM with 
the equilibrium liquid/gas transition in the canonical en¬ 
semble. Others include metastability, hysteresis, and the 
existence of critical nuclei. More generally, this anal¬ 
ogy suggests that the flocking transition should be seen 
as a phase-separation transition rather than an order- 
disorder transition. The fact that the liquid phase is 
ordered however plays a major role by forbidding a su¬ 
percritical region, which explains the atypical shape of 
the phase diagram. 

To construct a continuous theory for our model we first 
noticed that one needs to go beyond a standard mean- 
field approach. The latter indeed fails to capture the 
phase separation behavior because it lacks a density de¬ 
pendence of the transition temperature. Retaining part 
of the fluctuations neglected at the mean-field level then 
allowed us to derive a refined-mean-field model which 


accounts for the behavior of the microscopic model qual¬ 
itatively for all parameter values. 

The analytical solution for the phase-separated profile 
that we derived in sec.|V|is only one of a two-parameter 
family of solutions, as shown in m- Although it is the 
sole propagating solution accounting for phase separa¬ 
tion, the mecanism by which it is selected remains to 
be investigated. This is particularly interesting since, as 
shown in most of the picture laid out for the AIM re¬ 
mains valid for the Vicsek model, apart from the shape of 
the bands in the phase separated region. The full phase 
separation of the AIM is then replaced by a micro-phase 
separation, something which cannot be explained at the 
hydrodynamic level and necessit explicit noise terms. 

Beyond the sole case of the AIM, we showed that our 
results are also valid off lattice. We can thus consider the 
AIM as a representative example of a flocking model with 
discrete rotational symmetry. Variants with alignment 
between nearest neighbours, and not simply on-site, also 
yield similar results. 

Our study of the AIM relies on numerical simulations, 
microscopic derivation and study of hydrodynamic equa¬ 
tions. It says little about the universality of the emerging 
properties of the Active Ising Model and we strongly be¬ 
lieve that developing proper field theoretical approaches 
of the AIM and more general active spin models could 
shed light on a number of interesting questions. For in¬ 
stance, is the e = 0 limit of the AIM in the universality 
class of model C [50], which couples a conserved diffu¬ 
sive field and a non-conserved theory? Then, can one 
study the divergence of the correlation length of the AIM 
when approaching the T = 1 and v = 0 critical points? 
What are the corresponding universality classes? These 
questions will be addressed in future works. 

Last, the analogy of the phase transition in the AIM 
with an equilibrium liquid/gas transition triggers new 
questions. For example, could we define a mapping, at 
some ievei, with an equiiibrium system? And wouid it be 
possibie to change ensembie in this non-equiiibrium sys¬ 
tem, for exampie desining a grand-canonicai ensembie? 
These questions, if answered, wouid certainiy improve 
our theoreticai understanding of active matter systems. 
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Appendix A: One step beyond mean-field 

As shown in section |IV A[ the mean-field equations, 
which neglect all fluctuations and correlations, fail to de- 
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scribe the active Ising model since they predict a contin¬ 
uous phase transition between homogeneous phases. In 
this appendix we show how one can improve the mean- 
field approximation. To do so, we take into account 
the fluctuations of the local magnetisations and densities 
when computing the dynamics of their first moments (m) 
and (p). 


1. Gaussian fluctuations 

The simplest assumption that can be made about the 
fluctuations of m{x) and p{x) around their mean values, 
{m{x)) and (p(x)), is that they are Gaussian. For m(x), 
this can be seen as resulting from a central limit theo¬ 
rem: In a first approximation, the magnetisation is the 
sum of many spins fluctuating independently and, indeed. 
Fig. shows its fluctuations to be well described by a 
Gaussian. On the contrary, the distribution of the local 
density is not perfectly Gaussian, as shown in Fig. 

A better approximation could be obtained by consider¬ 
ing a Poisson distribution but, as will be apparent in the 
following, the first correction to mean-field comes from 
the fluctuations of m so this would not improve our ap¬ 
proximation. Furthermore we believe that, to improve 
our refined mean-field model, the next step should be to 
include the correlations between p and m, that we ne¬ 
glect in the following, and not higher cumulants of the 
distributions of p and m. 

More formally, the probability to observe a magneti¬ 
sation m and a density p at time t and position x given 
initial profiles po{x) and mo(x) are assumed to be given 
by 

V{p,m,x,t\po,mo) = N{p-p,(jl)N{m-fh,CTl^) (Al) 

where = e“®^/'^^/V 27 r(T 2 is the normal distri- 

bution and p(x, t) and fh{x^t) are the average value of 
the density and magnetisation fields. 

We further assume that the variances of the Gaus¬ 
sian distributions scales linearly with the local density: 

~ and = app{x,t). Again, the underly¬ 

ing assumption is that the fluctuations of the fields p{x) 
and m{x) arise from the sum of p independent contribu¬ 
tions. As shown in Fig. [^this is a rather good approx¬ 
imation in the gas phase, close to the critical point at 
/3 = I, p = oo. 


2. Corrections to mean-field 


In deriving hydrodynamic equations from Eq. (18) and 


(19), the only terms that have to be approximated are the 
non-linear contributions of the aligning interactions: 


l3m 


/3m \ 


/ = ( 2p sinh - -2m cosh - 

\ P P / 



^2/c+l \ 

I 



FIG. 23. Rescaled probability distributions of local density 
(left) and magnetisation (right) in the liquid phase at = I.I 
for different densities. ^(0,1) is the Gaussian distribution 
with zero mean and unit variance. D = I, £ = 0.9, L = 100. 



FIG. 24. Variance of the distribution of local density (left) 
and magnetisation (right) in the gas phase compared to a 
linear scaling. D = 1, £ = 0.9, L = 100. 


where 


Qj]^ — 2 


p2k-\-l 

{2k + l)\ 


p2k \ 

m) 


(A3) 


Using the assumption ( |A1| ), we can compute J as a 
sum of Gaussian integrals which can all be evaluated by 
saddle-point approximation in the limit of large p. We 
first notice that, since we neglect the correlations between 
p and m. 


(^) = {““+■> {Vp“> 


(A4) 


To compute we first change variables to = 

m — m so that 


/ -|-oo 

dm m^^+^A7(m — m, amp) 

-oo 

) 

du {u 3-m)‘^^~^^N'{u, amP) (A5) 


/ 


—oo 
+ 00 


We then expand in powers of u and compute the corre¬ 
sponding Gaussian integrals 


/ +00 2/c+l 

du ^ 

i=0 
k 


W(u, UmP) 




(A2) 


(A6) 
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Let us now evaluate the terms (p First, the inte¬ 
gral 

/*+oo 

= / dpp~‘^’^M{p-p,app) (A7) 

Jo 

is divergent because of the p = 0 lower limit. This is a 
simple discretisation problem which can be bypassed by 
introducing a cut-off ( at small density. For large p, the 
integrals will be dominated by large values of p so this 
cut-off does not play any role in the following. Changing 
variable to 5 = (p — p) /p, we find 

ds (1 + (A 8 ) 

y/27rap J4_i 

This integral can now be approximated by an asymptotic 
saddle-point expansion. In the limit of large p, the in¬ 
tegral is dominated by s 0. The lower limit of the 

integral ^ — 1 ^ —1 can thus be extended to — oo harm¬ 

lessly and one can expand ( 1 + 5 )“^^ to get the asymptotic 
expansion 


where 


a = -ai= p(l - 

(A16) 

_ aAi.iCo,! _ 3aam 

2 ~ 2 

(A17) 

3^ 

r2 = 3p{l3 - 3)amap + ^(/3 - 5)q:^ 

(A18) 


In practice, we take r 2 = 0 in the RMFM since the 
first order correction r/p suffices to a ccoun t for the phe¬ 
nomenology of the AIM. Expanding (AI4) to higher or¬ 
ders is not sufficient to get a quantitative agreement be¬ 
tween microscopic simulations and our refined mean-field 
model, probably because the most important correction 


to (AI5) would involve correlations between m and p. 
As we show in section however this first correction 
to mean-field is sufficient to capture the physics of the 
model. 


Appendix B: Linear stability analysis 




a^-2/c 




na 


2N 

E 

p i=0 

0{p 


‘Ik T i — I 


/ + ^ . - s2 

ds{—sye 

-oo 


-2N-2k- 


1 / 2 ) 


(A9) 

All the odd contributions vanish by symmetry. Changing 
variable to cj = ps‘^/{2ap), one recognises the integral 
form of a F function and finally 


N 


(^-“> = E 

i=o 


2k + 2j - 1\ r(j + I) 

U ) 




+ 0{p 


-2k-2N-l 


Putting everything together, we obtain 

oo k N 


(AlO) 


I — 


jy-^^-\-2k — 2i 


k=0i=0j=0 


) 2 /c+i 




I 


l+Ar+2/c-i > 


where 


bi.k = 


Cj,k 


2A: + I\ 2T(i-h 1/2) • 

2i ) 

2k + 2j - A 2^T{j + 1/2) 


2 j 




A 


ai 


(All) 

(A 12 ) 

(A13) 


Keeping only the dominant terms and reordering the sum 
in increasing powers of m yields 


^ = E 


^ ^l+2n p ^ / n h r 1 

i=0 j=0 


n=0 


p 2 n 


Expanding up to and I /p^, we finally obtain 

n3 


(A14) 


/~2(/3-l-1-(A15) 


The mean-field and refined mean-field equations read 
p = DAp — vdxTn (BI) 

Q 

Tfl 

rh = DAm — vdxp + 2mp — (B2) 

p^ 

where p = /3 — I — r/p and r = 0 for the mean-field equa¬ 
tions. These equations admit three steady homogeneous 
solutions p(x,t) = po, m(x,t) = mo- A disordered solu¬ 
tion with mo = 0 that exists for all po and /3, and two 
ordered solutions 


mo = ±po 



(B3) 


that exist only when p > 0 . 


1. Stability of the disordered profile 


Let us consider a small perturbation around the disor¬ 
dered profile, m(r, t) = ^m(r, t), p(r, t) = po -h^p(r, t). 
Going into Eourier space. 


6p 





and linearizing Eqs. 


and (B2), one finds 


f) {i 

*\6mJ \-iqxV -D\q\'^ + 2poJ \6mJ 


(B4) 


(B5) 


where we noted po = {j3 — 1 — r/po). The eigenvalues of 
the 2 x 2 matrix are 


A± 


“f ql) Mo 


Mo - 


(B 6 ) 
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FIG. 25. Real part of the largest eigenvalue A+ re lated to 
the stability of the ordered profile m = mo. of Eq. ( |B8[ ) for 
= r = = l and Qy = 0. Ordered profiles exist for 

all po > (fg = 2 but are unstable for (fg<po< (fi (red, green 
and yellow curves) and stable only for po > pn (blue curve). 
For the parameters considered here pa — 2.598. 


The profile is linearly unstable if one of these eigenval¬ 
ues has a positive real part. Clearly, the sign of /io 
controls the stability: the disordered profile is unsta¬ 
ble to long wavelength perturbations when /3 > 1 and 
PQ > = ^/(/^~ 1) and stable otherwise. This gives the 


first spinodal line pg in Fig. 16 


2. Stability of the ordered profile 


p(i', t) = po + (5/9(r, t) gives 


and the eigenvalues now read 


-iqxV \ f 5p\ 
-D\q\‘^ -4po) \SmJ 
(B7) 


A± = -D{ql+ql)-2pQ±. Upl - v'^ql 


2imoqxv{r + 2poPo) 


pI 


(B8) 

Equation ( |B8| ) shows qy to have a purely stabilizing ef¬ 
fect; Taking Qy = 0 thus does not affect the conclusions 
about the stability of the system. Computing numeri¬ 
cally 5R(A±) we observe (Fig. 25) that for small but pos¬ 
itive /io, 5R(A±) > 0 at long wave-length. The value of 
/io at which the system becomes stable can be deter¬ 
mined analytically as the point where dq^^{X±){qx = 0) 
changes sign (the first derivative being zero at = 0). 
This yields the second spinodal line shown in Fig. [T^ 


Pe = Pg 


vy^a {v‘^K + SD{(3 — 1)^) -h v‘^K. + SDa{l3 — 1) 


2v‘^i^ + SDa{P — 1 ) 


(B9) 


Linearizing the dynamics of a small perturbation 
around the ordered profile m(r, t) = mo + ^m(r,t), where k = 2a — 2(3. 
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